/*
  Copyright (C) 2011, Kyungjoo Kim
  
  This file is part of LINAL (LINear ALgebra)
  
  All rights reserved.

  This library is free software; you can redistribute it and/or
  modify it under the terms of the GNU Lesser General Public
  License as published by the Free Software Foundation; either
  version 2.1 of the License, or (at your option) any later version.
  
  This library is distributed in the hope that it will be useful,
  but WITHOUT ANY WARRANTY; without even the implied warranty of
  MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the GNU
  Lesser General Public License for more details.
  
  You should have received a copy of the GNU Lesser General Public
  License along with this library; if not, write to the Free Software
  Foundation, Inc., 51 Franklin Street, Fifth Floor, Boston, MA  02110-1301  USA
  Also add information on how to contact you by electronic and paper mail.
  
  
  Kyungjoo Kim
  iamkyungjoo@gmail.com
*/
#include "environment.hxx"

int main(int argc, char **argv) {

  linal::initialize(argc, argv);

  bool help;
  int datatype, nthreads, m, n;
  double tol;

  linal::get_option_bool  ("-help", "list options",
                           __FILE__,LINAL_OPTION_MODE_RUN, help);

  linal::get_option_int   ("-datatype", LINAL_OPTION_TEXT_DATATYPE,
                           __FILE__,TEST_DATATYPE, datatype, help);

  linal::get_option_int   ("-M", "size of matrix",
                           __FILE__,N, m, help);
  linal::get_option_int   ("-N", "size of matrix",
                           __FILE__,N, n, help);

  linal::get_option_int   ("-nthreads", "# of threads",
                           __FILE__,1, nthreads, help);

  linal::get_option_double("-tol", "error tolerence",
                           __FILE__,LINAL_ERROR_TOL, tol, help);

  if (help == LINAL_OPTION_MODE_DRY) {
    linal::finalize();
    return 0;
  }

  omp_set_num_threads(nthreads);

  linal::Flat_ A, B, d, D, norm;

  A.create( datatype, m, n );
  B.create( datatype, m, n );
  d.create( datatype, m, 1 );
  D.create( datatype, m, n );

  norm.create(linal::get_datatype_real(datatype), 1,1);

  FLA_Random_matrix( ~A );
  FLA_Copy( ~A, ~B );

  if (linal::get_verbose()) {
    FLA_Obj_show("= A =", ~A, "%4.3f", "===");
    FLA_Obj_show("= B =", ~B, "%4.3f", "===");
    FLA_Obj_show("= d =", ~d, "%4.3f", "===");
    FLA_Obj_show("= D =", ~D, "%4.3f", "===");
  }

  linal::host::internal::get_diag( ~A, ~d );
  linal::host::internal::set_diag( ~d, ~D );
  FLA_Set_diag( FLA_ZERO, ~A );

  if (linal::get_verbose()) {
    FLA_Obj_show("= A =", ~A, "%4.3f", "===");
    FLA_Obj_show("= B =", ~B, "%4.3f", "===");
    FLA_Obj_show("= d =", ~d, "%4.3f", "===");
    FLA_Obj_show("= D =", ~D, "%4.3f", "===");
  }

  linal::host::internal::get_diag( ~D, ~d );
  linal::host::internal::set_diag( ~d, ~A );

  FLA_Axpy( FLA_MINUS_ONE, ~A, ~B );
  FLA_Norm1( ~B, ~norm);

  if (linal::get_verbose()) 
    FLA_Obj_show("=Diff B=", ~B, "%4.3f", "===");

  float  norm_f = *((float *)FLA_Obj_buffer_at_view(~norm));
  double norm_d = *((double*)FLA_Obj_buffer_at_view(~norm));

  printf("- TEST::");
  for (int i=0;i<argc;++i)
    printf(" %s ", argv[i]);
  printf("\n");

  int rval;
  switch (linal::get_datatype_real(datatype)) {
  case LINAL_SINGLE_REAL:
    printf("SINGLE PRECISION\n");
    if (norm_f < LINAL_ERROR_TOL) {
      printf("PASS::Norm :: %10.9f \n", norm_f);   rval = 0;
    } else {
      printf("FAIL::Norm :: %10.9f \n", norm_f);   rval = -1;
    }
    break;
  case LINAL_DOUBLE_REAL:
    printf("DOUBLE PRECISION\n");
    if (norm_d < LINAL_ERROR_TOL) {
      printf("PASS::Norm :: %E \n", norm_d);   rval = 0;
    } else {
      printf("FAIL::Norm :: %E \n", norm_d);   rval = -1;
    }
    break;
  }

  // ---------------------------------------
  // ** Matrices  
  A.free();
  B.free();
  d.free();
  D.free();

  norm.free();

  linal::finalize();

  return rval;
}
